Environment-dependent metabolic investments in the mixotrophic chrysophyte Ochromonas

Code setup

Loading packages

Loads packages relevant to the script

# LOAD PACKAGES
require(Rmisc) # Has summarySE function

## Loading required package: Rmisc

## Loading required package: lattice

## Loading required package: plyr

require(sciplot) # Has bargraph.CI function

## Loading required package: sciplot

require(dismo) # For plotting convex hulls

## Loading required package: dismo

## Loading required package: raster

## Loading required package: sp

require(rworldmap) # For plotting map of locations

## Loading required package: rworldmap

## ### Welcome to rworldmap ###

## For a short introduction type :   vignette('rworldmap')

require(ape) # package for phylogenetic trees and Mantel test

## Loading required package: ape

## 
## Attaching package: 'ape'

## The following objects are masked from 'package:raster':
## 
##     rotate, zoom

Importing data

Imports data and performs some cleanup

# LOAD DATA
# NOTE: You will need to change the file path to wherever this spreadsheet is located. See instructions for a PC here: https://www.wikihow.com/Find-a-File%27s-Path-on-Windows

## PHYSIOLOGY
PAdata <- read.csv("/Users/hollyvm/GoogleSync/StudentWork/GinaBarbaglia/OchTradeoffData.csv")
#### PAdata <- read.csv("/Users/hollyvm/GoogleSync/StudentWork/GinaBarbaglia/PAdata_old.csv")


# Subset the data to exclude strain 2616 & light levels > 200
PAdata <- PAdata[PAdata$Strain!=2616 & PAdata$Light < 200 ,]


## GENETIC DISTANCE
gendist <- read.csv("/Users/hollyvm/GoogleSync/StudentWork/GinaBarbaglia/NumberOfDifferentBases.csv")


## STRAIN GEOGRAPHIC ORIGINS
metadata <- read.csv("/Users/hollyvm/GoogleSync/StudentWork/RyanMarczak/Ochromonas_TradeoffFunctions - CultureOrigins.csv")
metadata <- metadata[metadata$Strain!=2616,]

## BACTERIAL ABUNDANCE
bact.counts <- read.csv("/Users/hollyvm/GoogleSync/StudentWork/GinaBarbaglia/Ochromonas_TradeoffFunctions - PhotoAccBacterialCounts.csv")


## CN data
CNdata <- read.csv("/Users/hollyvm/GoogleSync/StudentWork/GinaBarbaglia/Ochromonas_TradeoffFunctions - Gina_CN.csv",header=T)

# PHOTOSYNTHESIS CALCULATIONS

Cfix.mult <- 0.01728178 # The multiplier 0.01728178 converts from electrons per chlorophyll molecule per second into pg C per pg chlorophyll-a per hour

PAdata$P_I.perC <- PAdata$chl.to.C*PAdata$P_I*Cfix.mult*12 # multiply by 12 because 12hrs light in each day; in mg of C per g C per day (because chl to C is in mg chl per g C)
PAdata$P_I.percell <- PAdata$chl.percell*PAdata$P_I*Cfix.mult*12 # multiply by 12 because 12hrs light in each day; in pg of C per cell per day (because chl per cell is in pg of chl per cell)

# Manually set photosynthetic rates to 0 if no light.
for(i in 1:dim(PAdata)[1]){
  if(PAdata$Light[i]==0){
    PAdata$P_I[i] <- 0
    PAdata$P_I.perC[i] <- 0
    PAdata$P_I.percell[i] <- 0 
    #### PAdata$P_I.cell[i] <- 0
  }
}

# GRAZING CALCULATIONS

fgC.per.bact <- 113.5996283 # from Lepori-Bui et al. 2022

# pg C of bacteria per pg C of Ochromonas per day = bacteria per och per day * 113 fgC/bact * 0.001 pg/fg * 1 och/pg C per och
PAdata$C.from.grazing.perCperday <- PAdata$Best.g*fgC.per.bact*0.001/PAdata$C.perOch

Setting plotting instructions and properties

# Rainbow of colours for the eight strains
# Colours from: https://github.com/BlakeRMills/MetBrewer - Hokusai 1 + a bit of Hokusai 2
col.a <- rgb(96/255,26/255,16/255)
col.b <- rgb(183/255,51/255,48/255)
col.c <- rgb(234/255,95/255,75/255)
col.d <- rgb(233/255,120/255,47/255)
col.e <- rgb(243/255,186/255,81/255)
col.f <- rgb(120/255,171/255,125/255)
col.g <- rgb(0/255,93/255,148/255)
col.h <- rgb(0/255,56/255,98/255)

# Rainbow of colours with reduced opacity for backgrounds of plots
alpha.bg <- .5
col.a.bg <- rgb(96/255,26/255,16/255,alpha.bg)
col.b.bg <- rgb(183/255,51/255,48/255,alpha.bg)
col.c.bg <- rgb(234/255,95/255,75/255,alpha.bg)
col.d.bg <- rgb(233/255,120/255,47/255,alpha.bg)
col.e.bg <- rgb(243/255,186/255,81/255,alpha.bg)
col.f.bg <- rgb(120/255,171/255,125/255,alpha.bg)
col.g.bg <- rgb(0/255,93/255,148/255,alpha.bg)
col.h.bg <- rgb(0/255,56/255,98/255,alpha.bg)

# Assembling into a vector of colors
col.strains <- c(col.a,col.b,col.c,col.d,col.e,col.f,col.g,col.h)
col.strains.bg <- c(col.a.bg,col.b.bg,col.c.bg,col.d.bg,col.e.bg,col.f.bg,col.g.bg,col.h.bg)

# "Mega" vector alternates saturated and unsaturated, for making barplots that show high and low bacterial treatments side by side
mega.col.strains <- rep(NaN,length(col.strains)*2)
for(i in 1:length(col.strains)){
  mega.col.strains[2*(i-1)+1] <- col.strains[i]
  mega.col.strains[2*i] <- col.strains.bg[i]
}

# Settings for symbol type and expansion factor of low (.X) and high (.R) bacterial treatments
pch.X <- 24
pch.R <- 21
ptcex.X <- .8
ptcex.R <- 1

Growth rates

Ochromonas growth rates increase with light and prey.

We have ordered this plot from left to right in terms of the strains that showed the biggest change in growth with the addition of bacteria. This chunk of code computes that.

# Figuring out strain ordering

# Make a dataframe summarizing the growth rate data
growth.summ <- summarySE(data=PAdata,'Growth.Rate',groupvars=c('Strain','Light','Bact.short'),na.rm=T)

# Compute the difference in growth rates between high and low bacterial treatments
growth.summ.x <- growth.summ[growth.summ$Bact.short=='X',]
growth.summ.r <- growth.summ[growth.summ$Bact.short=='R',]
growth.summ.x$growth.r <- growth.summ.r$Growth.Rate
growth.summ.x$delta.growth <- growth.summ.x$growth.r-growth.summ.x$Growth.Rate

# Make an additional summary across all light levels for each strain
uber.summ <- summarySE(data=growth.summ.x,'delta.growth',groupvars='Strain',na.rm=T)

# Sort the strains by the magnitude of this difference
uber.summ <- uber.summ[order(uber.summ$delta.growth,decreasing=T),]
strains <- uber.summ$Strain # the order command above enforces the order.

# Create a "Strain2" data column in PAdata that contains the strains as ordered; We'll use this later to make barplots in the same order.
PAdata$Strain2 <- factor(PAdata$Strain,levels=strains) 

This code creates a “rainbow plot function”, which we’ll use throughout the script

rainbowplot <- function(dataset,yvar,ylabel,xlabel,spacer,strains,linelevel,legendlocn,logchoice,ylimits,xlab.cex,toplabomit,xtickomit){
  
  # 1. Make the main plot
  if(missing(ylimits)){ # if no y-limits are imposed
  if(logchoice == 'y'){
    plot(dataset$Light,dataset[,colnames(dataset)==yvar],las=1,xlab='',ylab=ylabel, type = 'n',xlim=c(0,(150+spacer)*length(strains)-spacer),xaxt='n',log='y', ylim=c(min(dataset[,colnames(dataset)==yvar],na.rm=T),max(dataset[,colnames(dataset)==yvar],na.rm=T)))
  } else {
    plot(dataset$Light,dataset[,colnames(dataset)==yvar],las=1,xlab='',ylab=ylabel, type = 'n',xlim=c(0,(150+spacer)*length(strains)-spacer),xaxt='n') }
  } else { # if y-limits are imposed
  if(logchoice == 'y'){
    plot(dataset$Light,dataset[,colnames(dataset)==yvar],las=1,xlab='',ylab=ylabel, type = 'n',xlim=c(0,(150+spacer)*length(strains)-spacer),xaxt='n',log='y', ylim=ylimits)
  } else {
    plot(dataset$Light,dataset[,colnames(dataset)==yvar],las=1,xlab='',ylab=ylabel, type = 'n',xlim=c(0,(150+spacer)*length(strains)-spacer),xaxt='n',ylim=ylimits) }
    
  }
  
  # 2. Add an x-axis on the top (side 3) that indicates strain
  if(missing(toplabomit)){
  axis(side = 3, at = c(0:7)*(150+spacer)+75, labels=strains, lwd.ticks=0,line=-1,tck=0,lwd=0)
  }
  
  # 3. Label the x-axis
  if(missing(xlab.cex)){
  mtext(xlabel,side=1,line=1.5) # add x-axis label
  } else{
    mtext(xlabel,side=1,line=1.5,cex=xlab.cex)
  }

  # 4. Use a for loop to add a tick along the bottom for each strain
  for(i in 1:length(strains)){
    axis(side = 1, at = (i-1)*(150+spacer)+sort(unique(dataset$Light)),tck=.02,labels=F) # Add x-axis tick marks along the top that indicate light level
    if(missing(xtickomit)){
    axis(side=1, at = (i-1)*(150+spacer), tick=0, labels='0',line=-.75,cex.axis=1)
    axis(side=1, at = (i-1)*(150+spacer)+150, tick=0, labels='150',line=-.75,cex.axis=1)
  }}
  
  # 5. Add a horizontal dashed line to mark 0 growth
  abline(h=linelevel, lty=2)  
  
  # 6. Create a data summary
  dat.summ <- summarySE(data = dataset, measurevar = yvar, groupvars = c('Strain','Bact.short','Light'),na.rm=T)
  colnames(dat.summ)[5] <- 'YVar'
  
  # 7. Plot the data for each strain
for(i in 1:length(strains)){
  # 7.1. Subset the data
  subdat <- dataset[dataset$Strain==strains[i],]
  sub.dat.summ <- dat.summ[dat.summ$Strain==strains[i],]
  
  xcoord.adjust <- (i-1)*(150+spacer) # adjust x-coordinate to space out results
  
  # 7.2. Plot error bars for each light level
  arrows(sub.dat.summ$Light + xcoord.adjust, sub.dat.summ$YVar+sub.dat.summ$sd, sub.dat.summ$Light + xcoord.adjust, sub.dat.summ$YVar-sub.dat.summ$sd, length = 0.05, angle = 90, code= 3)
  
  # 7.3. Add lines connecting the points
  lines(sub.dat.summ[sub.dat.summ$Bact.short=='R',]$Light + xcoord.adjust,sub.dat.summ[sub.dat.summ$Bact.short=='R',]$YVar, lwd=2 , col = col.strains[i])
  lines(sub.dat.summ[sub.dat.summ$Bact.short=='X',]$Light + xcoord.adjust,sub.dat.summ[sub.dat.summ$Bact.short=='X',]$YVar, lwd=2 , col = col.strains[i])
  
  # 7.4. Add big points representing the mean values 
  points(sub.dat.summ[sub.dat.summ$Bact.short=='R',]$Light + xcoord.adjust,sub.dat.summ[sub.dat.summ$Bact.short=='R',]$YVar,pch=pch.R,cex=1.5*ptcex.R,bg=col.strains[i],col=col.strains[i])
  points(sub.dat.summ[sub.dat.summ$Bact.short=='X',]$Light + xcoord.adjust,sub.dat.summ[sub.dat.summ$Bact.short=='X',]$YVar,pch=pch.X,cex=1.5*ptcex.X,bg='white',col=col.strains[i])
}

  # 8. Add a legend
  legend(x = legendlocn, inset = 0.00, legend=c('High Bacteria', 'Low Bacteria'), pch=c(pch.R,pch.X), pt.bg=c('black','white'), horiz=T, pt.cex=1.5*c(ptcex.R,ptcex.X),cex=1,bty='n')

}

Call the function to make the plot of growth rates

# Set the physical amount of space between each strain's plot
spacer <- 80

# Use the 'rainbowplot' function defined above to make the plot. Syntax: 
#   First argument = dataset
#   Second argument = y-axis variable
#   Third argument = y-axis label
#   Fourth argument = x-axis label
#   Fifth argument = amount of space between plots
#   Sixth argument = order of the strains
#   Seventh argument = horizontal line position, or NA if none needed
#   Eigth argument = location for the legend
#   Ninth argument = log scaling of y axis, enter 'y' or 'n'
#   Tenth argument = y axis limits (optional)
#   Eleventh argument = expansion factor for x-axis label (optional)
#   Twelfth argument = 'n' if omitting top labels (optional)
#   xtickomit == omits x tick labels (optional)
rainbowplot(PAdata,'Growth.Rate','Growth Rate (per day)','Light Level (grouped by strain)',spacer,strains,0,'topright','n')

Adding significance tests for positive growth rates at zero light.

for(i in 1:length(strains)){
StrainChoice <- strains[i]

print(StrainChoice)
print(t.test(PAdata[PAdata$Strain==StrainChoice & PAdata$Light==0 & PAdata$Bact.short.2=='X',]$Growth.Rate,alternative='greater'))
print(t.test(PAdata[PAdata$Strain==StrainChoice & PAdata$Light==0 & PAdata$Bact.short.2=='X',]$Growth.Rate,alternative='greater')$p.value*8)

print(t.test(PAdata[PAdata$Strain==StrainChoice & PAdata$Light==0 & PAdata$Bact.short.2=='R',]$Growth.Rate,alternative='greater'))
print(t.test(PAdata[PAdata$Strain==StrainChoice & PAdata$Light==0 & PAdata$Bact.short.2=='R',]$Growth.Rate,alternative='greater')$p.value*8)

}

## [1] 584
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "X", ]$Growth.Rate
## t = -23.044, df = 2, p-value = 0.9991
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.0700703        Inf
## sample estimates:
## mean of x 
##  -0.06219 
## 
## [1] 7.992489
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = 11.114, df = 2, p-value = 0.004
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  0.2226941       Inf
## sample estimates:
## mean of x 
## 0.3020567 
## 
## [1] 0.03199761
## [1] 590
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "X", ]$Growth.Rate
## t = -2.1186, df = 2, p-value = 0.9159
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.02265696         Inf
## sample estimates:
##    mean of x 
## -0.009526667 
## 
## [1] 7.326883
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = 28.202, df = 2, p-value = 0.0006275
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  0.2482359       Inf
## sample estimates:
## mean of x 
## 0.2769067 
## 
## [1] 0.005019862
## [1] 2951
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "X", ]$Growth.Rate
## t = -19.717, df = 2, p-value = 0.9987
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.04668147         Inf
## sample estimates:
## mean of x 
##  -0.04066 
## 
## [1] 7.989751
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = 29.145, df = 2, p-value = 0.0005876
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  0.2075958       Inf
## sample estimates:
## mean of x 
##   0.23071 
## 
## [1] 0.00470064
## [1] 1393
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "X", ]$Growth.Rate
## t = 14.334, df = 2, p-value = 0.002416
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  0.01637981        Inf
## sample estimates:
## mean of x 
##   0.02057 
## 
## [1] 0.01932594
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = 21.994, df = 2, p-value = 0.00103
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  0.09395351        Inf
## sample estimates:
## mean of x 
## 0.1083367 
## 
## [1] 0.008243506
## [1] 1392
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "X", ]$Growth.Rate
## t = -1.9748, df = 2, p-value = 0.9065
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.0944444        Inf
## sample estimates:
##   mean of x 
## -0.03810333 
## 
## [1] 7.252083
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = 1.0883, df = 2, p-value = 0.1951
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.02722462         Inf
## sample estimates:
##  mean of x 
## 0.01617667 
## 
## [1] 1.560467
## [1] 1148
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "X", ]$Growth.Rate
## t = -3.2879, df = 2, p-value = 0.9593
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.05593171         Inf
## sample estimates:
##   mean of x 
## -0.02962333 
## 
## [1] 7.67451
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = 11.701, df = 2, p-value = 0.003613
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  0.1581005       Inf
## sample estimates:
## mean of x 
## 0.2106767 
## 
## [1] 0.02890121
## [1] 1150
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "X", ]$Growth.Rate
## t = -7.0857, df = 2, p-value = 0.9903
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.0858412        Inf
## sample estimates:
## mean of x 
##  -0.06079 
## 
## [1] 7.922634
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = 4.4165, df = 2, p-value = 0.02382
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  0.004724663         Inf
## sample estimates:
##  mean of x 
## 0.01394333 
## 
## [1] 0.1905366
## [1] 1391
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "X", ]$Growth.Rate
## t = -21.356, df = 2, p-value = 0.9989
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.05259256         Inf
## sample estimates:
##   mean of x 
## -0.04626667 
## 
## [1] 7.991259
## 
##  One Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == 0 & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = 1.1713, df = 2, p-value = 0.1811
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  -0.03581775         Inf
## sample estimates:
## mean of x 
##   0.02399 
## 
## [1] 1.448596

Significance tests for Strain 1391 growth rates

lights <- unique(PAdata$Light)

StrainChoice <- 1391

for(i in 1:length(lights)){

print(lights[i])
print(t.test(PAdata[PAdata$Strain==StrainChoice & PAdata$Light==lights[i] & PAdata$Bact.short.2=='X',]$Growth.Rate,PAdata[PAdata$Strain==StrainChoice & PAdata$Light==lights[i] & PAdata$Bact.short.2=='R',]$Growth.Rate,alternative='two.sided'))
print(t.test(PAdata[PAdata$Strain==StrainChoice & PAdata$Light==lights[i] & PAdata$Bact.short.2=='X',]$Growth.Rate,PAdata[PAdata$Strain==StrainChoice & PAdata$Light==lights[i] & PAdata$Bact.short.2=='R',]$Growth.Rate,alternative='two.sided')$p.value*7)

}

## [1] 25
## 
##  Welch Two Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "X", ]$Growth.Rate and PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = -0.10613, df = 2.0278, p-value = 0.925
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1763517  0.1677517
## sample estimates:
## mean of x mean of y 
## 0.4416667 0.4459667 
## 
## [1] 6.475311
## [1] 75
## 
##  Welch Two Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "X", ]$Growth.Rate and PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = -1.8876, df = 2.4393, p-value = 0.1763
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2098331  0.0664998
## sample estimates:
## mean of x mean of y 
## 0.5282333 0.5999000 
## 
## [1] 1.234052
## [1] 150
## 
##  Welch Two Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "X", ]$Growth.Rate and PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = -3.6938, df = 2.6564, p-value = 0.04216
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.259405458 -0.009661209
## sample estimates:
## mean of x mean of y 
## 0.5352667 0.6698000 
## 
## [1] 0.2950884
## [1] 50
## 
##  Welch Two Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "X", ]$Growth.Rate and PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = -7.2012, df = 2.8018, p-value = 0.006908
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.11583895 -0.04282772
## sample estimates:
## mean of x mean of y 
## 0.4692667 0.5486000 
## 
## [1] 0.04835458
## [1] 100
## 
##  Welch Two Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "X", ]$Growth.Rate and PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = -0.50281, df = 2.1783, p-value = 0.6614
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1911777  0.1483110
## sample estimates:
## mean of x mean of y 
## 0.5826333 0.6040667 
## 
## [1] 4.629555
## [1] 15
## 
##  Welch Two Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "X", ]$Growth.Rate and PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = -0.99345, df = 3.4645, p-value = 0.3847
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.07523450  0.03736783
## sample estimates:
## mean of x mean of y 
## 0.3571333 0.3760667 
## 
## [1] 2.692982
## [1] 0
## 
##  Welch Two Sample t-test
## 
## data:  PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "X", ]$Growth.Rate and PAdata[PAdata$Strain == StrainChoice & PAdata$Light == lights[i] & PAdata$Bact.short.2 == "R", ]$Growth.Rate
## t = -3.4111, df = 2.0447, p-value = 0.07392
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.15704345  0.01653012
## sample estimates:
##   mean of x   mean of y 
## -0.04626667  0.02399000 
## 
## [1] 0.5174078

# Set the physical amount of space between each strain's plot
spacer <- 80

# Use the 'rainbowplot' function defined above to make the plot. Syntax: 
#   First argument = dataset
#   Second argument = y-axis variable
#   Third argument = y-axis label
#   Fourth argument = x-axis label
#   Fifth argument = amount of space between plots
#   Sixth argument = order of the strains
#   Seventh argument = horizontal line position, or NA if none needed
#   Eigth argument = location for the legend
#   Ninth argument = log scaling of y axis, enter 'y' or 'n'
#   Tenth argument = y axis limits (optional)
#   Eleventh argument = expansion factor for x-axis label (optional)
#   Twelfth argument = 'n' if omitting top labels (optional)
#   xtickomit == omits x tick labels (optional)
rainbowplot(PAdata,'Growth.Rate','Growth Rate (per day)','Light Level (grouped by strain)',spacer,strains,0,'topright','n')

horiz.adj <- 30
x.coords.left <- (c(1:length(strains))-1)*(150+spacer)-horiz.adj
x.coords.right <- x.coords.left + 2*horiz.adj

# Add significance values for rice treatments
signif.rice <- c('*','**','**','**','','*','','')
text(x.coords.left,rep(.005,length(strains)),signif.rice,cex=1.5,srt=90,pos=4)

# Add significance values for xenic treatments
signif.xenic <- c('','','','*','','','','')
text(x.coords.right,rep(.005,length(strains)),signif.xenic,cex=1.5,col='gray50',pos=4,srt=90)

Growth increases and saturates with increasing photosynthetic rates, and is generally higher with higher bacterial availability.

Computing carbon use efficiency

Carbon use efficiency (CUE) = amount of carbon used for growth / (total amount of carbon acquired)

We have the growth rate, which is in units of per day. Or, more explicitly, in units of carbon turned into new biomass per carbon of current biomass per day. That’s the numerator of CUE.

For the carbon acquired, we need to sum up the carbon gained through photosynthesis, and the carbon gained through grazing.

For photosynthesis, this is photosynthetic rate * amount of photosynthetic machinery * 12 hrs of sunlight per day. But we have to be careful to do all of this in units of per gram of carbon.

# Growth.Rate <- rate of growth, in units of per day
# C.perOch <- ng C per Ochromonas cell

# chl.percell <- pg chl-a per cell
# chl.to.C <- mg chl-a per g C
# P_I <- electrons transferred per molecule of chlorophyll-a 

Cfix.mult <- 0.01728178 # The multiplier 0.01728178 converts from electrons per chlorophyll molecule per second into pg C per pg chlorophyll-a per hour

PAdata$chl.to.C <- PAdata$chl.percell/PAdata$C.perOch*1000 # mg chl-a per g carbon

# pg C fixed via photosynthesis per pg of Ochromonas per day = electrons per molecule of chl-a per second * [pg C per pg chl-a per hr / electrons per molecule of chl-a per second] * 12hrs of light per day * pg chl per pg C
PAdata$P_I.perC <- PAdata$P_I*Cfix.mult*12*PAdata$chl.to.C/1000 # multiply by 12 because 12hrs light in each day; divide by 1000 because chl to c in mg per g


# Best.g <- number of bacteria estimated to be eaten per Ochromonas cell per day

fgC.per.bact <- 113.5996283 # fentograms C per bacterial cell
# pg C of bacteria per pg C of Ochromonas per day = bacteria per och per day * 113 fgC/bact * 0.001 pg/fg * 1 och/pg C per och
PAdata$C.from.grazing.perCperday <- PAdata$Best.g*fgC.per.bact*0.001/PAdata$C.perOch # scaled to g bacterial C per g Ochromonas C per day

PAdata$attack.per.C.scaled <- PAdata$attack.per.C*24*1000*1000 # scaled to mL per mg carbon per day

PAdata$CUE <- PAdata$Growth.Rate/(PAdata$C.from.grazing.perCperday+PAdata$P_I.perC)

Make a barplot of the results

# Make the barplot
bargraph.CI(Strain2,CUE,group=Bact.short, data=PAdata,las=1,legend=F,xlab='Ochromonas Strain',ylab='Carbon Use Efficiency',leg.lab=c('High Bacteria','Low Bacteria'),col=mega.col.strains,err.width = .05,ylim=c(0,1.199))
legend(x='topleft',inset=0.0,legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=2.3,pch=22,pt.bg=c('black','gray80'),horiz=T)

# Uncomment to see results of the Tukey test
#TukeyHSD(aov(CUE~Bact.short*Strain2,data=PAdata))

# Add text labels to indicate significance
text(x = c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5,16.5,17.5,19.5,20.5,22.5,23.5),y=c(0.35,0.85,0.43,.7,.44,.72,.325,.64,.32,.815,.85,.925,.69,1.18,.25,.635),c('a','b','a','b','a','b','a','b','a','b','b','b','b','c','a','b'))

# Make the barplot
bargraph.CI(Strain2,CUE,group=Bact.short, data=PAdata,las=1,legend=F,xlab='Ochromonas Strain',ylab='Carbon Use Efficiency',leg.lab=c('High Bacteria','Low Bacteria'),col=c('gray50','white'),err.width = .05,ylim=c(0,1.199))
legend(x='topleft',inset=0.0,legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=2.3,pch=22,pt.bg=c('gray50','white'),horiz=T)

# Uncomment to see results of the Tukey test
#TukeyHSD(aov(CUE~Bact.short*Strain2,data=PAdata))

# Add text labels to indicate significance
text(x = c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5,16.5,17.5,19.5,20.5,22.5,23.5),y=c(0.35,0.85,0.43,.7,.44,.72,.325,.64,.32,.815,.85,.925,.69,1.18,.25,.635),c('a','b','a','b','a','b','a','b','a','b','b','b','b','c','a','b'))

# Make the barplot
mega.col.strains2 <- rep(NaN,length(col.strains)*2)
for(i in 1:length(col.strains)){
  mega.col.strains2[2*(i-1)+1] <- col.strains[i]
  mega.col.strains2[2*i] <- 'white'
}

mega.col.strains3 <- rep(NaN,length(col.strains)*2)
for(i in 1:length(col.strains)){
  mega.col.strains3[2*(i-1)+1] <- col.strains[i]
  mega.col.strains3[2*i] <-  col.strains[i]
}
bargraph.CI(Strain2,CUE,group=Bact.short, data=PAdata,las=1,legend=F,xlab='Ochromonas Strain',ylab='Carbon Use Efficiency',leg.lab=c('High Bacteria','Low Bacteria'),col=mega.col.strains2,border='black',err.width = .05,ylim=c(0,1.199))
legend(x='topleft',inset=c(0.01,0.01),legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=2.3,pch=22,pt.bg=c('gray20','white'),horiz=T)
#bargraph.CI(Strain2,CUE,group=Bact.short, data=PAdata,las=1,legend=F,err.width = .05,add=T)

# Uncomment to see results of the Tukey test
#TukeyHSD(aov(CUE~Bact.short*Strain2,data=PAdata))

# Add text labels to indicate significance
text(x = c(1.5,2.5,4.5,5.5,7.5,8.5,10.5,11.5,13.5,14.5,16.5,17.5,19.5,20.5,22.5,23.5),y=c(0.35,0.85,0.43,.7,.44,.72,.325,.64,.32,.815,.85,.925,.69,1.18,.25,.635),c('a','b','a','b','a','b','a','b','a','b','b','b','b','c','a','b'))

We were curious about whether carbon use efficiency tended to decline as growth rates increased. Here’s a plot of CUE by strain as a function of growth rate. Generally, no support for this hypothesis.

par(mar=c(2,2.5,1.5,1),mfrow=c(2,4))
for(i in 1:length(strains)){
  subdat <- PAdata[PAdata$Strain==strains[i]&PAdata$Light<200,]
  plot(subdat$Growth.Rate,subdat$CUE,xlab='Growth Rate',ylab='Carbon Use Efficiency',main=strains[i],las=1,pch=c(21,22)[as.factor(subdat$Bact.short)],cex=1.5,col=col.strains[i],bg=c(col.strains.bg[i],'white')[as.factor(subdat$Bact.short)])
}

Explaining growth rates with photosynthesis and phagotrophy

par(mar=c(4,4,1,1),mfrow=c(1,1))

# PHOTOSYNTHESIS PREDICTS GROWTH: PHOTOSYNTHESIS PER CARBON

plot(PAdata$Growth.Rate~PAdata$P_I.perC,col=col.strains[as.factor(PAdata$Strain)],las=1,xlab='Photosynthetic Rate (g C per g C per d)', ylab='Growth Rate (per day)',lwd=.5,pch=21,type='n')

mu.summ <- summarySE(data = PAdata, measurevar = 'Growth.Rate', groupvars = c('Strain','Bact.short','Light'),na.rm=T)
PI.summ <- summarySE(data = PAdata, measurevar = 'P_I.perC', groupvars = c('Strain','Bact.short','Light'),na.rm=T)

for(i in 1:length(strains)){
  sub.mu <- mu.summ[mu.summ$Strain==strains[i],]
  sub.pi <- PI.summ[PI.summ$Strain==strains[i],]
  arrows(sub.pi$P_I.perC,sub.mu$Growth.Rate+sub.mu$se,sub.pi$P_I.perC,sub.mu$Growth.Rate-sub.mu$se,code=3,angle=90,length=0.05)
  arrows(sub.pi$P_I.perC+sub.pi$se,sub.mu$Growth.Rate,sub.pi$P_I.perC-sub.pi$se,sub.mu$Growth.Rate,code=3,angle=90,length=0.05)
  points(sub.mu$Growth.Rate~sub.pi$P_I.perC,pch=c(pch.R,pch.X)[as.factor(sub.mu$Bact.short)],col=col.strains[i],lwd=2,cex=1.5*c(ptcex.R,ptcex.X)[as.factor(sub.mu$Bact.short)],bg=c(col.strains[i],'white')[as.factor(sub.mu$Bact.short)])
}

rect(-10,-10,100,10,col=rgb(1,1,1,.5)) # Adds a white rectangle that's 50% transparent to make the points "fade out" a little bit

for(i in 1:length(strains)){
  sub.mu <- mu.summ[mu.summ$Strain==strains[i],]
  sub.pi <- PI.summ[PI.summ$Strain==strains[i],]
  
  data.x <- as.data.frame(cbind(sub.mu[sub.mu$Bact.short=='X',]$Growth.Rate,sub.pi[sub.pi$Bact.short=='X',]$P_I.perC))
  colnames(data.x) <- c('Growth','PI')
  c.est <- data.x[1,1]
  a.est <- 5
  data.x[1,2] <- 0
  fit.x <- nls(Growth~a*PI/(h + PI)+c, data = data.x, start = list(a = a.est, h =1, c = c.est))
  pi.min <- min(data.x$PI,na.rm=T); pi.max <- max(data.x$PI,na.rm=T); pi.set <- seq(from=pi.min,to = pi.max,length.out=50)
  mu.set <- summary(fit.x)$parameters[1,1]*pi.set / (summary(fit.x)$parameters[2,1]+pi.set)+summary(fit.x)$parameters[3,1]
  lines(pi.set,mu.set,lwd=2,col=col.strains[i],lty=2)
  
  data.r <- as.data.frame(cbind(sub.mu[sub.mu$Bact.short=='R',]$Growth.Rate,sub.pi[sub.pi$Bact.short=='R',]$P_I.perC))
  colnames(data.r) <- c('Growth','PI')
  c.est <- data.r[1,1]
  a.est <- (data.r[2,1]-data.r[1,1])/data.r[2,2]
  fit.r <- nls(Growth~a*PI/(h + PI)+c, data = data.r, start = list(a = a.est, h =1, c = c.est))
  pi.min <- min(data.r$PI,na.rm=T); pi.max <- max(data.r$PI,na.rm=T); pi.set <- seq(from=pi.min,to = pi.max,length.out=50)
  mu.set <- summary(fit.r)$parameters[1,1]*pi.set / (summary(fit.r)$parameters[2,1]+pi.set)+summary(fit.r)$parameters[3,1]
  lines(pi.set,mu.set,lwd=2,col=col.strains[i])
}

legend(x = 'bottomright', inset = 0.01, legend=strains, pch=21,pt.cex=1.5, col=col.strains,pt.bg=col.strains, title='Strain',pt.lwd=2)
legend(x = .7, y = 0.1, legend=c('High','Low'), title='Prey',pch=c(pch.R,pch.X),pt.cex=1.5*c(ptcex.R,ptcex.X),pt.bg=c('black','white'),pt.lwd=2)

Revision of this plot dividing by strain

spacer <- 80/150

 # 1. Make the main plot
plot(1,1,las=1,xlab='',ylab='Growth Rate (per day)', type = 'n',xlim=c(0,(1+spacer)*length(strains)-spacer),xaxt='n',ylim=c(-.1,1))
  
# 2. Add an x-axis on the top (side 3) that indicates strain
axis(side = 3, at = c(0:7)*(1+spacer)+0.5, labels=strains, lwd.ticks=0,line=-1,tck=0,lwd=0)
  
# 3. Label the x-axis
mtext('Photosynthetic Rate (g C per g C per day)',side=1,line=1.5) 

# 4. Use a for loop to add a tick along the bottom for each strain
for(i in 1:length(strains)){
  axis(side = 1, at = (i-1)*(1+spacer)+c(0,0.25,0.5,0.75,1),tck=.02,labels=F) # Add x-axis tick marks along the top that indicate light level
  axis(side=1, at = (i-1)*(1+spacer), tick=0, labels='0',line=-.75,cex.axis=1)
  axis(side=1, at = (i-1)*(1+spacer)+1, tick=0, labels='1',line=-.75,cex.axis=1)
  axis(side=1, at = (i-1)*(1+spacer)+.5, tick=0, labels='0.5',line=-.75,cex.axis=1)
  }
  
# 5. Add a horizontal dashed line to mark 0 growth
abline(h=0, lty=2)  

# 6. Summarize the data
mu.summ <- summarySE(data = PAdata, measurevar = 'Growth.Rate', groupvars = c('Strain','Bact.short','Light'),na.rm=T)
PI.summ <- summarySE(data = PAdata, measurevar = 'P_I.perC', groupvars = c('Strain','Bact.short','Light'),na.rm=T)
  
# 7. Plot the data for each strain
for(i in 1:length(strains)){
  # 7.1. Adjust x-axis plotting
  xcoord.adjust <- (i-1)*(1+spacer) # adjust x-coordinate to space out results
  
  # 7.2. Subset the data
  sub.mu <- mu.summ[mu.summ$Strain==strains[i],]
  sub.pi <- PI.summ[PI.summ$Strain==strains[i],]
  
  # 7.3. Make xenic curve fits
  # Subset and give recognizable column names
  data.x <- as.data.frame(cbind(sub.mu[sub.mu$Bact.short=='X',]$Growth.Rate,sub.pi[sub.pi$Bact.short=='X',]$P_I.perC)); colnames(data.x) <- c('Growth','PI')
  # Estimate parameters
  c.est <- data.x[1,1]
  a.est <- 5
  data.x[1,2] <- 0 # Double-check that photosynthetic rate in darkness = 0
  # Perform curve fit
  fit.x <- nls(Growth~a*PI/(h + PI)+c, data = data.x, start = list(a = a.est, h =1, c = c.est))
  # Generate predictions
  pi.min <- min(data.x$PI,na.rm=T); pi.max <- max(data.x$PI,na.rm=T); pi.set <- seq(from=pi.min,to = pi.max,length.out=50)
  mu.set <- summary(fit.x)$parameters[1,1]*pi.set / (summary(fit.x)$parameters[2,1]+pi.set)+summary(fit.x)$parameters[3,1]
  # Plot a line for the prediction
  lines(pi.set+xcoord.adjust,mu.set,lwd=2,col=col.strains[i],lty=2)
  
  # 7.4. Make rice curve fits
  data.r <- as.data.frame(cbind(sub.mu[sub.mu$Bact.short=='R',]$Growth.Rate,sub.pi[sub.pi$Bact.short=='R',]$P_I.perC))
  colnames(data.r) <- c('Growth','PI')
  c.est <- data.r[1,1]
  a.est <- (data.r[2,1]-data.r[1,1])/data.r[2,2]
  fit.r <- nls(Growth~a*PI/(h + PI)+c, data = data.r, start = list(a = a.est, h =1, c = c.est))
  pi.min <- min(data.r$PI,na.rm=T); pi.max <- max(data.r$PI,na.rm=T); pi.set <- seq(from=pi.min,to = pi.max,length.out=50)
  mu.set <- summary(fit.r)$parameters[1,1]*pi.set / (summary(fit.r)$parameters[2,1]+pi.set)+summary(fit.r)$parameters[3,1]
  lines(pi.set+xcoord.adjust,mu.set,lwd=2,col=col.strains[i])
  
  # 7.5. Plot x-error bars
  arrows(sub.pi$P_I.perC+sub.pi$se+xcoord.adjust,sub.mu$Growth.Rate,sub.pi$P_I.perC-sub.pi$se+xcoord.adjust,sub.mu$Growth.Rate,code=3,angle=90,length=0.05)
  # 7.6. Plot y-error bars
  arrows(sub.pi$P_I.perC+xcoord.adjust,sub.mu$Growth.Rate+sub.mu$se,sub.pi$P_I.perC+xcoord.adjust,sub.mu$Growth.Rate-sub.mu$se,code=3,angle=90,length=0.05)
  # 7.7. Plot big points representing the mean values
  points(sub.pi$P_I.perC+xcoord.adjust,sub.mu$Growth.Rate,pch=c(pch.R,pch.X)[as.factor(sub.mu$Bact.short)],col=col.strains[i],lwd=2,cex=1.5*c(ptcex.R,ptcex.X)[as.factor(sub.mu$Bact.short)],bg=c(col.strains[i],'white')[as.factor(sub.mu$Bact.short)])
  
}

  # 8. Add a legend
  legend(x = 'topright', inset = 0.00, legend=c('High Bacteria', 'Low Bacteria'), pch=c(pch.R,pch.X), pt.bg=c('black','white'), horiz=T, pt.cex=1.5*c(ptcex.R,ptcex.X),cex=1,bty='n')

Let’s try it with 8 distinct panels so that we can expand the x-axis and see the saturating shape.

par(mar=c(1,1,0.1,0.1))
layout(matrix(c(9,1,1,1,2,2,2,3,3,3,4,4,4,9,1,1,1,2,2,2,3,3,3,4,4,4,9,1,1,1,2,2,2,3,3,3,4,4,4,9,5,5,5,6,6,6,7,7,7,8,8,8,9,5,5,5,6,6,6,7,7,7,8,8,8,9,5,5,5,6,6,6,7,7,7,8,8,8,rep(10,13)),nrow=7,ncol=13,byrow = T))

subplotlabels <- c("(a)","(b)",'(c)',"(d)","(e)","(f)","(g)","(h)")

# Curve fit data storage

phot.fit <- as.data.frame(rep(strains,each = 2))
colnames(phot.fit) <- 'strains'
phot.fit$bact <- rep(c('X','R'),length(strains))
phot.fit$c <- NaN
phot.fit$a <- NaN
phot.fit$a.err <- NaN


# 6. Summarize the data
mu.summ <- summarySE(data = PAdata, measurevar = 'Growth.Rate', groupvars = c('Strain','Bact.short','Light'),na.rm=T)
PI.summ <- summarySE(data = PAdata, measurevar = 'P_I.perC', groupvars = c('Strain','Bact.short','Light'),na.rm=T)
  
# 7. Plot the data for each strain
for(i in 1:length(strains)){
  plot(1,1,type='n',xlab='Phot',ylab='Growth',las=1,ylim=c(-.050,1),xlim=c(0,1.1),xaxt='n',yaxt='n')
  axis(side=1,at=c(0,.2,.4,.6,.8,1),labels=F)
  axis(side=2,at=c(0,.2,.4,.6,.8,1),labels=F)
  if(i %in% c(1,5)){axis(side=2,at=c(0,.2,.4,.6,.8,1),labels=T,las=1,cex.axis=1.2)}
  if(i > 4){axis(side=1,at=c(0,.2,.4,.6,.8,1),labels=T,las=1,cex.axis=1.2)}
  
  # 7.2. Subset the data
  sub.mu <- mu.summ[mu.summ$Strain==strains[i],]
  sub.pi <- PI.summ[PI.summ$Strain==strains[i],]
  
  # 7.3. Make xenic curve fits
  # Subset and give recognizable column names
  data.x <- as.data.frame(cbind(sub.mu[sub.mu$Bact.short=='X',]$Growth.Rate,sub.pi[sub.pi$Bact.short=='X',]$P_I.perC)); colnames(data.x) <- c('Growth','PI')
  # Estimate parameters
  c.est <- data.x[1,1]
  a.est <- 5
  data.x[1,2] <- 0 # Double-check that photosynthetic rate in darkness = 0
  # Perform curve fit
  fit.x <- nls(Growth~a*PI/(h + PI)+c, data = data.x, start = list(a = a.est, h =1, c = c.est))
  # Generate predictions
  pi.min <- min(data.x$PI,na.rm=T); pi.max <- max(data.x$PI,na.rm=T); pi.set <- seq(from=pi.min,to = pi.max,length.out=50)
  mu.set <- summary(fit.x)$parameters[1,1]*pi.set / (summary(fit.x)$parameters[2,1]+pi.set)+summary(fit.x)$parameters[3,1]
  # Plot a line for the prediction
  lines(pi.set,mu.set,lwd=2,col=col.strains[i],lty=2)
  
  phot.fit$c[2*i-1] <- summary(fit.x)$parameters[3,1]
  phot.fit$a[2*i-1] <- summary(fit.x)$parameters[1,1]
  phot.fit$a.err[2*i-1] <- summary(fit.x)$parameters[1,2]
  
  # 7.4. Make rice curve fits
  data.r <- as.data.frame(cbind(sub.mu[sub.mu$Bact.short=='R',]$Growth.Rate,sub.pi[sub.pi$Bact.short=='R',]$P_I.perC))
  colnames(data.r) <- c('Growth','PI')
  c.est <- data.r[1,1]
  a.est <- (data.r[2,1]-data.r[1,1])/data.r[2,2]
  fit.r <- nls(Growth~a*PI/(h + PI)+c, data = data.r, start = list(a = a.est, h =1, c = c.est))
  pi.min <- min(data.r$PI,na.rm=T); pi.max <- max(data.r$PI,na.rm=T); pi.set <- seq(from=pi.min,to = pi.max,length.out=50)
  mu.set <- summary(fit.r)$parameters[1,1]*pi.set / (summary(fit.r)$parameters[2,1]+pi.set)+summary(fit.r)$parameters[3,1]
  lines(pi.set,mu.set,lwd=2,col=col.strains[i])
  
  #points(1,summary(fit.r)$parameters[1,1]+summary(fit.r)$parameters[3,1])
  
  phot.fit$c[2*i] <- summary(fit.r)$parameters[3,1]
  phot.fit$a[2*i] <- summary(fit.r)$parameters[1,1]
  phot.fit$a.err[2*i] <- summary(fit.r)$parameters[1,2]
  
  # 7.5. Plot x-error bars
  arrows(sub.pi$P_I.perC+sub.pi$se,sub.mu$Growth.Rate,sub.pi$P_I.perC-sub.pi$se,sub.mu$Growth.Rate,code=3,angle=90,length=0.05)
  # 7.6. Plot y-error bars
  arrows(sub.pi$P_I.perC,sub.mu$Growth.Rate+sub.mu$se,sub.pi$P_I.perC,sub.mu$Growth.Rate-sub.mu$se,code=3,angle=90,length=0.05)
  # 7.7. Plot big points representing the mean values
  points(sub.pi$P_I.perC,sub.mu$Growth.Rate,pch=c(pch.R,pch.X)[as.factor(sub.mu$Bact.short)],col=col.strains[i],lwd=2,cex=1.5*c(ptcex.R,ptcex.X)[as.factor(sub.mu$Bact.short)],bg=c(col.strains[i],'white')[as.factor(sub.mu$Bact.short)])
  
  # Add strain label
  text(1.15,0,paste('CCMP ',strains[i],sep=''),cex=1.5,col=col.strains[i],pos=2)
  
  # Add panel label
  text(0.05,.95,subplotlabels[i],cex=1.5)
  
  # Add legend
  if(i == 4){ legend(x = .71, y = 1.01, legend=c('High','Low'), title='Prey',pch=c(pch.R,pch.X),pt.cex=1.5*c(ptcex.R,ptcex.X),pt.bg=c('black','white'),pt.lwd=2,cex=1.2)
 }
  
}

# Add overall axis labels
text(x = -1.4, y = -.33, "Photosynthetic Rate (g C per g C per day)",srt=0,xpd=NA,cex=1.5)
text(x = -4.2, y = 1.15, "Growth Rate (per day)",srt=90,xpd=NA,cex=1.5)

Supplemental figure showing saturations

fill.col <- rep('white',16)
for(i in 1:length(strains)){
  fill.col[2*i-1] <- col.strains.bg[i]
  fill.col[2*i] <- col.strains[i]
}

phot.fit$a.alt <- phot.fit$a
phot.fit$a.err.alt <- phot.fit$a.err
phot.fit[phot.fit$strains==2951 & phot.fit$bact=='R',]$a.alt <- phot.fit[phot.fit$strains==2951 & phot.fit$bact=='R',]$a-12.5
#phot.fit[phot.fit$strains==2951 & phot.fit$bact=='R',]$a.err.alt <- .2

barplot(phot.fit$a.alt+phot.fit$c,las=1,ylim=c(0,2.8),space=c(.3,.1),border='black',col=fill.col,xlab='',ylab='',yaxt='n')
axis(side=2,at=c(0,.5,1,1.5,2,2.5),labels=c(0,.5,1,1.5,14.5,15),las=1)
mtext('Estimated Saturating Growth Rate (per day)',side=2,line=2.8)
mtext('Strain x Bacterial Treatment',side=1,line=1.3)

x.coords <- c(0.8,1.9,3.2,4.3,5.6,6.7,8,9.1,10.4,11.5,12.8,13.9,15.2,16.3,17.6,18.7)
arrows(x.coords,phot.fit$a.alt+phot.fit$c + phot.fit$a.err.alt,x.coords, phot.fit$a.alt + phot.fit$c - phot.fit$a.err.alt,code=3,angle=90,length=0.05)

axis(side=1, at = c(1.35,3.75,6.15,8.55,10.95,13.35,15.75,18.15),labels=strains,tick=F,line=-1,col=col.strains)

lines(c(-1,1),c(1.7,1.8)+.01,lwd=6,col='white',xpd=T)
#lines(c(-1,1),c(1.7,1.8)+.0,lwd=1,col='black',xpd=T)
lines(c(-.7,-.3)+.03,c(1.7,1.73)-.005,lwd=1,col='black',xpd=T)
lines(c(-.7,-.3)+.03,c(1.7,1.73)+.045,lwd=1,col='black',xpd=T)
lines(c(6,8),c(1.7,1.7)+.05,lwd=8,col=rgb(1,1,1,.5),xpd=T)
lines(c(6,8),c(1.7,1.7)+.05,lwd=6,col=rgb(1,1,1,.5),xpd=T)
lines(c(6,8),c(1.7,1.7)+.05,lwd=5,col=rgb(1,1,1,1),xpd=T)
lines(c(6,8),c(1.7,1.7)+.05,lwd=4,col=rgb(1,1,1,1),xpd=T)

text(c(1.35,3.75,6.15,8.55,10.95,13.35,15.75,18.15),y=c(1.2,1,2.4,1.3,0,0,0,1),c('*','*','','*','','','','*'),cex=1.3)

legend(x='topright',inset=c(0.1,0.1),legend=c('Low Bacteria','High Bacteria'),pch=22,pt.cex=1.5,pt.bg=c('gray80','gray20'))

Supplemental figure for heterotrophy

par(mar=c(4,4,1,1),mfrow=c(1,1))

plot(PAdata$Growth.Rate~PAdata$C.from.grazing.perCperday,col=col.strains[as.factor(PAdata$Strain)],las=1,xlab='Phagotrophy Rate (g C per g C per d)', ylab='Growth Rate (per day)',lwd=.5,pch=21,type='n',log='x',xlim=c(0.01,50),ylim=c(0.15,1))

mu.summ <- summarySE(data = PAdata, measurevar = 'Growth.Rate', groupvars = c('Strain','Bact.short','Light'),na.rm=T)
het.summ <- summarySE(data = PAdata, measurevar = 'C.from.grazing.perCperday', groupvars = c('Strain','Bact.short','Light'),na.rm=T)

for(i in 1:length(strains)){
  sub.mu <- mu.summ[mu.summ$Strain==strains[i],]
  sub.pi <- het.summ[het.summ$Strain==strains[i],]
  arrows(sub.pi$C.from.grazing.perCperday,sub.mu$Growth.Rate+sub.mu$se,sub.pi$C.from.grazing.perCperday,sub.mu$Growth.Rate-sub.mu$se,code=3,angle=90,length=0.05)
  arrows(sub.pi$C.from.grazing.perCperday+sub.pi$se,sub.mu$Growth.Rate,sub.pi$C.from.grazing.perCperday-sub.pi$se,sub.mu$Growth.Rate,code=3,angle=90,length=0.05)
  points(sub.mu$Growth.Rate~sub.pi$C.from.grazing.perCperday,pch=c(pch.R,pch.X)[as.factor(sub.mu$Bact.short)],col=col.strains[i],lwd=2,cex=1.5*c(ptcex.R,ptcex.X)[as.factor(sub.mu$Bact.short)],bg=c(col.strains[i],'white')[as.factor(sub.mu$Bact.short)])
}

legend(x = 'topleft', inset = 0.01, legend=strains, pch=21,pt.cex=1.5, col=col.strains,pt.bg=col.strains,title='Strain',pt.lwd=2)
#legend(x = 'topleft', inset = 0.01, legend=strains, pch=21,pt.cex=1.5, col=col.strains,pt.bg=col.strains,title='Strain')

legend(x = .04, y = 1.025, legend=c('High','Low'), title='Prey',pch=c(pch.R,pch.X),pt.cex=1.5*c(ptcex.R,ptcex.X),pt.bg=c('black','white'),pt.lwd=2)

Phenotypic plasticity

Mixotroph metabolic investment traits vary with environmental conditions. (1) Chlorophyll investment decreases with light. (2) Attack rates decrease with prey availability

# Set the amount of space between each curve
spacer <- 80

# Set up 2-row figure

matrixreps <- 10
layout(matrix(c(3,rep(1,matrixreps),rep(2,matrixreps),4),nrow=2+2*matrixreps,ncol=1))
par(mar=c(1,4,0,0.1))


# TOP ROW: CHLOROPHYLL
rainbowplot(PAdata,'chl.to.C','Chlorophyll (mg Chl-a per g C)','',spacer,strains,NA,'topright','n',xtickomit='T')
mtext('(a)',side=1,at=-235,line=-16)

# BOTTOM ROW: ATTACK RATE
rainbowplot(PAdata[PAdata$Light>0,],'attack.per.C.scaled','Attack rate (mL per mg C per day)','Light Level (grouped by strain)',spacer,strains,NA,'topleft','y',c(.08,20),toplabomit='T',xlab.cex = .7)
mtext('(b)',side=1,at=-235,line=-16)

Tradeoffs

Some lineages experience sharp tradeoffs between investments: As chlorophyll per carbon increases, attack rates per carbon decrease. These tradeoffs only seem to emerge when growth rates are low, implying that cells are resource-limited and thus might be maximizing their investments in each form of metabolism.

# Compile summary data

PAdata$Strain.Bact.Light <- paste(PAdata$Strain,PAdata$Bact.short,PAdata$Light,sep='.')
growth.summ.C <- summarySE(data=PAdata[PAdata$Light>0,],'Growth.Rate',groupvars=c('Strain.Bact.Light'),na.rm=T)
chl.summ.C <- summarySE(data=PAdata[PAdata$Light>0,],'chl.to.C',groupvars=c('Strain.Bact.Light'),na.rm=T)
att.summ.C <- summarySE(data=PAdata[PAdata$Light>0,],'attack.per.C.scaled',groupvars=c('Strain.Bact.Light'),na.rm=T)


growth.summ.C$mean.chl <- chl.summ.C$chl.to.C[match(growth.summ.C$Strain.Bact.Light,chl.summ.C$Strain.Bact.Light)]
growth.summ.C$se.chl <- chl.summ.C$se[match(growth.summ.C$Strain.Bact.Light,chl.summ.C$Strain.Bact.Light)]

growth.summ.C$mean.att <- att.summ.C$attack.per.C.scaled[match(growth.summ.C$Strain.Bact.Light,att.summ.C$Strain.Bact.Light)]
growth.summ.C$se.att <- att.summ.C$se[match(growth.summ.C$Strain.Bact.Light,att.summ.C$Strain.Bact.Light)]


growth.summ.C$Strain <- PAdata$Strain[match(growth.summ.C$Strain.Bact.Light,PAdata$Strain.Bact.Light)]
growth.summ.C$Bact.short <- PAdata$Bact.short[match(growth.summ.C$Strain.Bact.Light,PAdata$Strain.Bact.Light)]


# Set up plot
par(mar=c(2,2.2,.4,.4),mfrow=c(2,4))
subplotlabels <- c("(a)","(b)",'(c)',"(d)","(e)","(f)","(g)","(h)")

for(i in 1:length(strains)){
  plot(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$attack.per.C.scaled,xlab='Photo. Invest.', ylab='Hetero. Invest.',type='n',las=1,cex.axis=1.5)
  
  # Subset the data for the strain
  subdat <- growth.summ.C[growth.summ.C$Strain==strains[i],]
  
  # Set up color ramp based on growth rate
  subdat$order <- findInterval(subdat$Growth.Rate,sort(subdat$Growth.Rate))
  pal = colorRampPalette(c("white",col.strains[i]))
  #pal = colorRampPalette(c(col.strains[i],"white"))
  
  # Add convex hull
  ch <- convHull(cbind(subdat$mean.chl,subdat$mean.att))
  plot(ch,add=T,col=pal(10)[2],lty=3,border='gray50')
  
  # Add error bars
  arrows(subdat$mean.chl+subdat$se.chl,subdat$mean.att,subdat$mean.chl-subdat$se.chl,subdat$mean.att,angle=90,code=3,length=.05)
  arrows(subdat$mean.chl,subdat$mean.att+subdat$se.att,subdat$mean.chl,subdat$mean.att-subdat$se.att,angle=90,code=3,length=.05)
  
  # Add points
  points(subdat$mean.chl,subdat$mean.att,pch=c(21,24)[as.factor(subdat$Bact.short)],cex=1.3*c(1.6,1.5)[as.factor(subdat$Bact.short)],lwd=1,bg=pal(nrow(subdat))[subdat$order]) #c('goldenrod','dodgerblue')[as.factor(sub.chl$Bact.short)])
  text(mean(c(min(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,na.rm=T),max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,na.rm=T))),max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$attack.per.C.scaled,na.rm=T)*.92,paste('CCMP ',strains[i],sep=''),col=col.strains[i],pos=3,cex=1.5)
  text(min(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,na.rm=T)*1.08,max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$attack.per.C.scaled,na.rm=T)*.92,subplotlabels[i],pos=3,cex=1.5)
}

Vertical version:

# Compile summary data

PAdata$Strain.Bact.Light <- paste(PAdata$Strain,PAdata$Bact.short,PAdata$Light,sep='.')
growth.summ.C <- summarySE(data=PAdata[PAdata$Light>0,],'Growth.Rate',groupvars=c('Strain.Bact.Light'),na.rm=T)
chl.summ.C <- summarySE(data=PAdata[PAdata$Light>0,],'chl.to.C',groupvars=c('Strain.Bact.Light'),na.rm=T)
att.summ.C <- summarySE(data=PAdata[PAdata$Light>0,],'attack.per.C.scaled',groupvars=c('Strain.Bact.Light'),na.rm=T)


growth.summ.C$mean.chl <- chl.summ.C$chl.to.C[match(growth.summ.C$Strain.Bact.Light,chl.summ.C$Strain.Bact.Light)]
growth.summ.C$se.chl <- chl.summ.C$se[match(growth.summ.C$Strain.Bact.Light,chl.summ.C$Strain.Bact.Light)]

growth.summ.C$mean.att <- att.summ.C$attack.per.C.scaled[match(growth.summ.C$Strain.Bact.Light,att.summ.C$Strain.Bact.Light)]
growth.summ.C$se.att <- att.summ.C$se[match(growth.summ.C$Strain.Bact.Light,att.summ.C$Strain.Bact.Light)]


growth.summ.C$Strain <- PAdata$Strain[match(growth.summ.C$Strain.Bact.Light,PAdata$Strain.Bact.Light)]
growth.summ.C$Bact.short <- PAdata$Bact.short[match(growth.summ.C$Strain.Bact.Light,PAdata$Strain.Bact.Light)]


# Set up plot
par(mar=c(2,2.2,.4,.4),mfrow=c(4,2))
subplotlabels <- c("(a)","(b)",'(c)',"(d)","(e)","(f)","(g)","(h)")

for(i in 1:length(strains)){
  plot(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$attack.per.C.scaled,xlab='Photo. Invest.', ylab='Hetero. Invest.',type='n',las=1,cex.axis=1.5)
  
  # Subset the data for the strain
  subdat <- growth.summ.C[growth.summ.C$Strain==strains[i],]
  
  # Set up color ramp based on growth rate
  subdat$order <- findInterval(subdat$Growth.Rate,sort(subdat$Growth.Rate))
  pal = colorRampPalette(c("white",col.strains[i]))
  #pal = colorRampPalette(c(col.strains[i],"white"))
  
  # Add convex hull
  ch <- convHull(cbind(subdat$mean.chl,subdat$mean.att))
  plot(ch,add=T,col=pal(10)[2],lty=3,border='gray50')
  
  # Add error bars
  arrows(subdat$mean.chl+subdat$se.chl,subdat$mean.att,subdat$mean.chl-subdat$se.chl,subdat$mean.att,angle=90,code=3,length=.05)
  arrows(subdat$mean.chl,subdat$mean.att+subdat$se.att,subdat$mean.chl,subdat$mean.att-subdat$se.att,angle=90,code=3,length=.05)
  
  # Add points
  points(subdat$mean.chl,subdat$mean.att,pch=c(21,24)[as.factor(subdat$Bact.short)],cex=1.3*c(1.6,1.5)[as.factor(subdat$Bact.short)],lwd=1,bg=pal(nrow(subdat))[subdat$order]) #c('goldenrod','dodgerblue')[as.factor(sub.chl$Bact.short)])
  y.coordinate <- max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$attack.per.C.scaled,na.rm=T)-(max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$attack.per.C.scaled,na.rm=T)-min(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$attack.per.C.scaled,na.rm=T))*.1
  text(mean(c(min(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,na.rm=T),max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,na.rm=T))),y.coordinate,paste('CCMP ',strains[i],sep=''),col=col.strains[i],pos=3,cex=1.5)
  text(min(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,na.rm=T)*1.08,y.coordinate,subplotlabels[i],pos=3,cex=1.5)
}

Partitioning by light

par(mar=c(4,4,1,1),mfrow=c(2,3))
LLs <- c(15,25,50,75,100,150)
growth.summ.C$Light <- PAdata$Light[match(growth.summ.C$Strain.Bact.Light,PAdata$Strain.Bact.Light)]

for(i in 1:length(LLs)){
  plot(PAdata[PAdata$Light==LLs[i],]$chl.to.C,PAdata[PAdata$Light==LLs[i],]$attack.per.C.scaled,xlab='Photo. Invest.', ylab='Hetero. Invest.',type='n',las=1,cex.axis=1.5)
  
  # Subset the data for the strain
  subdat <- growth.summ.C[growth.summ.C$Light==LLs[i],]
  
  # Set up color ramp based on growth rate
  subdat$order <- findInterval(subdat$Growth.Rate,sort(subdat$Growth.Rate))
  pal = colorRampPalette(c("white",col.strains[i]))
  #pal = colorRampPalette(c(col.strains[i],"white"))
  
  # Add convex hull
  #ch <- convHull(cbind(subdat$mean.chl,subdat$mean.att))
  #plot(ch,add=T,col=pal(10)[2],lty=3,border='gray50')
  
  # Add error bars
  arrows(subdat$mean.chl+subdat$se.chl,subdat$mean.att,subdat$mean.chl-subdat$se.chl,subdat$mean.att,angle=90,code=3,length=.05)
  arrows(subdat$mean.chl,subdat$mean.att+subdat$se.att,subdat$mean.chl,subdat$mean.att-subdat$se.att,angle=90,code=3,length=.05)
  

  # Add points

  #points(subdat$mean.chl,subdat$mean.att,pch=c(21,24)[as.factor(subdat$Bact.short)],cex=1.3*c(1.6,1.5)[as.factor(subdat$Bact.short)],lwd=1,bg=pal(nrow(subdat))[subdat$order]) #c('goldenrod','dodgerblue')[as.factor(sub.chl$Bact.short)])
  #y.coordinate <- max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$attack.per.C.scaled,na.rm=T)-(max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$attack.per.C.scaled,na.rm=T)-min(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$attack.per.C.scaled,na.rm=T))*.1
  #text(mean(c(min(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,na.rm=T),max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,na.rm=T))),y.coordinate,paste('CCMP ',strains[i],sep=''),col=col.strains[i],pos=3,cex=1.5)
  #text(min(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$chl.to.C,na.rm=T)*1.08,y.coordinate,subplotlabels[i],pos=3,cex=1.5)
  
    for(j in 1:length(strains)){
    subsubdat <- subdat[subdat$Strain==strains[j]&subdat$Bact.short=='X',]
    points(subsubdat$mean.chl,subsubdat$mean.att,pch=21,bg=col.strains[j],cex=2)
  }
}

Measuring synergies vs tradeoffs

Just because you have a tradeoff in how you can make investments doesn’t mean that you can’t experience synergies in the yields of those investments.

What would this look like?

We know how to look for tradeoffs in investments by looking for cases in which attack rates and chlorophyll have negative relationships. But how can we observe synergies?

1. Photosynthetic efficiency, rate of photosynthesis per chlorophyll, and rate of photosynthesis per carbon in high and low prey levels

# Set up 3-row figure
par(mar=c(4,4,1,1),mfrow=c(3,1))

layout(matrix(c(4,1,1,1,2,2,2,3,3,3,5),nrow=11,ncol=1))
par(mar=c(1,4,0,.3))

## ROW 1:

# Set the amount of space between each curve
spacer <- 80

# TOP ROW: Photosynthetic efficiency
#rainbowplot(PAdata[PAdata$Light>-10,],'Fv.Fm',expression('Photosynthetic Efficiency (Fv/Fm)'),'',spacer,strains,NA,'bottomleft','n',xlab.cex=.7,xtickomit = T)
#mtext('(a)',side=1,at=-180,line=-16.6)
rainbowplot(PAdata[PAdata$Light>-10,],'Fv.Fm',substitute(paste('Photosynthetic Efficiency (',italic('F'),'v/',italic('F'),"m)"))
,'',spacer,strains,NA,'bottomleft','n',xlab.cex=.7,xtickomit = T)
mtext('(a)',side=1,at=-180,line=-16.6)

# MIDDLE ROW: PHOTOSYNTHESIS PER CHLOROPHYLL
rainbowplot(PAdata,'P_I','Photosynthetic Rate (e- per chl-a per s)','',spacer,strains,NA,'topleft','n',xlab.cex=.7,xtickomit = T,toplabomit = T)
mtext('(b)',side=1,at=-180,line=-16.6)

# BOTTOM ROW: PHOTOSYNTHESIS PER CARBON
rainbowplot(PAdata,'P_I.perC','Photosynthetic Rate (g C per g C per d)','Light Level (grouped by strain)',spacer,strains,NA,'topright','n',xlab.cex=.7,toplabomit = T)
mtext('(c)',side=1,at=-180,line=-15.9)

2. Rates of grazing as a function of photosynthetic rates

# Add new summary data

phot.summ.C <- summarySE(data=PAdata[PAdata$Light>0,],'P_I.perC',groupvars=c('Strain.Bact.Light'),na.rm=T)
graz.summ.C <- summarySE(data=PAdata[PAdata$Light>0,],'C.from.grazing.perCperday',groupvars=c('Strain.Bact.Light'),na.rm=T)


growth.summ.C$mean.phot <- phot.summ.C$P_I.perC[match(growth.summ.C$Strain.Bact.Light,phot.summ.C$Strain.Bact.Light)]
growth.summ.C$se.phot <- phot.summ.C$se[match(growth.summ.C$Strain.Bact.Light,phot.summ.C$Strain.Bact.Light)]

growth.summ.C$mean.graz <- graz.summ.C$C.from.grazing.perCperday[match(growth.summ.C$Strain.Bact.Light,graz.summ.C$Strain.Bact.Light)]
growth.summ.C$se.graz <- graz.summ.C$se[match(growth.summ.C$Strain.Bact.Light,graz.summ.C$Strain.Bact.Light)]


# Set up plot
#par(mar=c(4,4,1,1),mfrow=c(2,4))
par(mar=c(2,2.7,.4,.3),mfrow=c(2,4))


for(i in 1:length(strains)){
  plot(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$P_I.perC,PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$C.from.grazing.perCperday,xlab='Photo. Rate', ylab='Hetero. Rate',type='n',las=1,cex.axis=1.5)
  
  # Subset the data for the strain
  subdat <- growth.summ.C[growth.summ.C$Strain==strains[i],]
  
  # Set up color ramp based on growth rate
  subdat$order <- findInterval(subdat$Growth.Rate,sort(subdat$Growth.Rate))
  pal = colorRampPalette(c("white",col.strains[i]))
  #pal = colorRampPalette(c(col.strains[i],"white"))
  
  # Add convex hull
  ch <- convHull(cbind(subdat$mean.phot,subdat$mean.graz))
  plot(ch,add=T,col=pal(10)[2],lty=3,border='gray50')
  
  # Add error bars
  arrows(subdat$mean.phot+subdat$se.phot,subdat$mean.graz,subdat$mean.phot-subdat$se.phot,subdat$mean.graz,angle=90,code=3,length=.05)
  arrows(subdat$mean.phot,subdat$mean.graz+subdat$se.graz,subdat$mean.phot,subdat$mean.graz-subdat$se.graz,angle=90,code=3,length=.05)
  
  # Add points
  points(subdat$mean.phot,subdat$mean.graz,pch=c(21,24)[as.factor(subdat$Bact.short)],cex=1.3*c(1.6,1.5)[as.factor(subdat$Bact.short)],lwd=1,bg=pal(nrow(subdat))[subdat$order])
  
  text(mean(c(min(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$P_I.perC,na.rm=T),max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$P_I.perC,na.rm=T))),max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$C.from.grazing.perCperday,na.rm=T)*.925,paste('CCMP ',strains[i],sep=''),col=col.strains[i],pos=3, cex=1.5)
  text(min(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$P_I.perC,na.rm=T)*1.15,max(PAdata[PAdata$Light>0&PAdata$Strain==strains[i],]$C.from.grazing.perCperday,na.rm=T)*.92,subplotlabels[i],pos=3,cex=1.5)
}

Do a nitrogen budget

First let’s make bargraphs of the C per cell, N per cell, and C:N ratios of the Ochromonas.

# Some data cleanup
CNdata$Strain2 <- factor(CNdata$Strain,levels=strains) # Orders the strains for plotting according to growth differential
CNdata$Bact.short <- substr(CNdata$Bact,1,1) # Uses the first character of the bacterial treatment as the identifier
CNdata <- CNdata[CNdata$Bact.short%in%c("X","R"),] # Eliminates additional data (e.g., CN analyzer controls and empty rows)

layout(matrix(c(1,1,1,2,2,2,3,3,3,4),nrow=10))
par(mar=c(1,5,2,.5))

bargraph.CI(Strain2,C.perOch,group=Bact.short, data=CNdata,las=1,legend=F,xlab='',ylab='',col=mega.col.strains,err.width = .05, cex.axis=1.5, cex.names=1.5)
legend(x='topleft',inset=0.0,legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=3,pch=22,pt.bg=c('gray30','gray80'),cex=1.5)
mtext('pg C per cell',side=2,line=3,cex=1.)
mtext('(a)',side=1,line=-16.5,at=-2,cex=1.3)

bargraph.CI(Strain2,N.perOch,group=Bact.short, data=CNdata,las=1,legend=F,xlab='',ylab='',col=mega.col.strains,err.width = .05, cex.axis=1.5, cex.names=1.5)
#legend(x='topleft',inset=0.0,legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=3,pch=22,pt.bg=c('gray30','gray80'),cex=1.5)
mtext('pg N per cell',side=2,line=3,cex=1.)
mtext('(b)',side=1,line=-16.5,at=-2,cex=1.3)

bargraph.CI(Strain2,CN.Och,group=Bact.short, data=CNdata,las=1,legend=F,xlab='Ochromonas Strain',ylab='',col=mega.col.strains,err.width = .05, cex.axis=1.5, cex.names=1.5)
#legend(x='topleft',inset=0.0,legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=3,pch=22,pt.bg=c('gray30','gray80'),cex=1.5)
mtext('Ochromonas Strain',side=1,line=3,cex=1)
mtext('C:N Ratio',side=2,line=3,cex=1.)
mtext('(c)',side=1,line=-16.5,at=-2,cex=1.3)

Significance testing

# Carbon per cell
for(i in 1:length(strains)){
  strainchoice <- strains[i]
  print(strainchoice)
  print(t.test(CNdata[CNdata$Strain2==strainchoice&CNdata$Bact.short=='X',]$C.perOch,CNdata[CNdata$Strain2==strainchoice&CNdata$Bact.short=='R',]$C.perOch)$p.value*8)
}

## [1] 584
## [1] 1.612919
## [1] 590
## [1] 0.1427455
## [1] 2951
## [1] 0.04330356
## [1] 1393
## [1] 1.131674
## [1] 1392
## [1] 2.885748
## [1] 1148
## [1] 0.004714536
## [1] 1150
## [1] 0.002260486
## [1] 1391
## [1] 0.2408857

# Nitrogen per cell
for(i in 1:length(strains)){
  strainchoice <- strains[i]
  print(strainchoice)
  print(t.test(CNdata[CNdata$Strain2==strainchoice&CNdata$Bact.short=='X',]$N.perOch,CNdata[CNdata$Strain2==strainchoice&CNdata$Bact.short=='R',]$N.perOch)$p.value*8)
}

## [1] 584
## [1] 0.9298087
## [1] 590
## [1] 0.03041344
## [1] 2951
## [1] 0.05840129
## [1] 1393
## [1] 0.2628777
## [1] 1392
## [1] 0.4966082
## [1] 1148
## [1] 0.006270415
## [1] 1150
## [1] 8.928302e-05
## [1] 1391
## [1] 0.1179823

# C:N ratio per cell
for(i in 1:length(strains)){
  strainchoice <- strains[i]
  print(strainchoice)
  print(t.test(CNdata[CNdata$Strain2==strainchoice&CNdata$Bact.short=='X',]$CN.Och,CNdata[CNdata$Strain2==strainchoice&CNdata$Bact.short=='R',]$CN.Och)$p.value*8)
}

## [1] 584
## [1] 0.8608805
## [1] 590
## [1] 0.3552752
## [1] 2951
## [1] 0.00310438
## [1] 1393
## [1] 0.04421683
## [1] 1392
## [1] 0.07194664
## [1] 1148
## [1] 0.01494959
## [1] 1150
## [1] 2.523798
## [1] 1391
## [1] 0.003594524

Adding points based on light level.

# Some data cleanup
CNdata$Strain2 <- factor(CNdata$Strain,levels=strains) # Orders the strains for plotting according to growth differential
CNdata$Bact.short <- substr(CNdata$Bact,1,1) # Uses the first character of the bacterial treatment as the identifier
CNdata <- CNdata[CNdata$Bact.short%in%c("X","R"),] # Eliminates additional data (e.g., CN analyzer controls and empty rows)

layout(matrix(c(1,1,1,2,2,2,3,3,3,4),nrow=10))
par(mar=c(1,5,2,.5))

bargraph.CI(Strain2,C.perOch,group=Bact.short, data=CNdata,las=1,legend=F,xlab='',ylab='',col=mega.col.strains,err.width = .05, cex.axis=1.5, cex.names=1.5,ylim=c(0,70))
legend(x='topleft',inset=0.0,legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=3,pch=22,pt.bg=c('gray30','gray80'),cex=1.5)
mtext('pg C per cell',side=2,line=3,cex=1.)
mtext('(a)',side=1,line=-16.5,at=-2,cex=1.3)

x.coords <- c(2,5,8,11,14,17,20,23)

for(i in 1:length(strains)){
  # Add xenic results
  subdat <- CNdata[CNdata$Strain==strains[i]&CNdata$Bact=='X',]
  points(rep(x.coords[i]+.5,dim(subdat)[1]),subdat$C.perOch,pch=c(21,24))
  # Add rice results
  subdat <- CNdata[CNdata$Strain==strains[i]&CNdata$Bact=='R',]
  points(rep(x.coords[i]-.5,dim(subdat)[1]),subdat$C.perOch,pch=c(21,24))
  
}

text(x.coords,y=c(0,0,43,0,0,67,44,0),c('','','*','','','**','**',''),cex=2,lwd=2.5)

bargraph.CI(Strain2,N.perOch,group=Bact.short, data=CNdata,las=1,legend=F,xlab='',ylab='',col=mega.col.strains,err.width = .05, cex.axis=1.5, cex.names=1.5, ylim=c(0,11))
#legend(x='topleft',inset=0.0,legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=3,pch=22,pt.bg=c('gray30','gray80'),cex=1.5)
mtext('pg N per cell',side=2,line=3,cex=1.)
mtext('(b)',side=1,line=-16.5,at=-2,cex=1.3)
for(i in 1:length(strains)){
  # Add xenic results
  subdat <- CNdata[CNdata$Strain==strains[i]&CNdata$Bact=='X',]
  points(rep(x.coords[i]+.5,dim(subdat)[1]),subdat$N.perOch,pch=c(21,24))
  # Add rice results
  subdat <- CNdata[CNdata$Strain==strains[i]&CNdata$Bact=='R',]
  points(rep(x.coords[i]-.5,dim(subdat)[1]),subdat$N.perOch,pch=c(21,24))
  
}
text(x.coords,y=c(0,7.8,7,0,0,9.5,6,0),c('','*','*','','','**','***',''),cex=2,lwd=2.5)

bargraph.CI(Strain2,CN.Och,group=Bact.short, data=CNdata,las=1,legend=F,xlab='Ochromonas Strain',ylab='',col=mega.col.strains,err.width = .05, cex.axis=1.5, cex.names=1.5, ylim=c(0,25))
#legend(x='topleft',inset=0.0,legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=3,pch=22,pt.bg=c('gray30','gray80'),cex=1.5)
mtext('Ochromonas Strain',side=1,line=3,cex=1)
mtext('C:N Ratio',side=2,line=3,cex=1.)
mtext('(c)',side=1,line=-16.5,at=-2,cex=1.3)
for(i in 1:length(strains)){
  # Add xenic results
  subdat <- CNdata[CNdata$Strain==strains[i]&CNdata$Bact=='X',]
  points(rep(x.coords[i]+.5,dim(subdat)[1]),subdat$CN.Och,pch=c(21,24))
  # Add rice results
  subdat <- CNdata[CNdata$Strain==strains[i]&CNdata$Bact=='R',]
  points(rep(x.coords[i]-.5,dim(subdat)[1]),subdat$CN.Och,pch=c(21,24))
  
}
text(x.coords,y=c(0,0,12,21,0,23,0,13),c('','','**','*','','*','','**'),cex=2,lwd=2.5)

Here we are calculating the amount of nitrogen that the cells get from grazing, vs how much they need to get in order to make their own machinery. Negative numbers mean they have to be taking up some inorganic N directly.

CN.bacteria <- 106/16 # 10.1 # Average number from Michelle's data. 
PAdata$N.from.grazing.perCperday <- PAdata$C.from.grazing.perCperday/CN.bacteria # Convert carbon from bacteria into nitrogen from bacteria
PAdata$N.from.grazing.perNperday <- PAdata$N.from.grazing.perCperday*(PAdata$C.perOch/PAdata$N.perOch) #multiply by Ochromonas C:N ratio to get Nitrogen per Ochromonas Nitrogen per day

PAdata$N.deficit <- (PAdata$N.from.grazing.perNperday - PAdata$Growth.Rate)/PAdata$Growth.Rate

bargraph.CI(Strain2,N.deficit,group=Bact.short, data=PAdata[PAdata$Light>0,],las=1,legend=F,xlab='Ochromonas Strain',ylab='N into biomass per N acq via grazing',leg.lab=c('High Bacteria','Low Bacteria'),col=mega.col.strains,err.width = .05)
legend(x='topleft',inset=0.05,legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=2,pch=22,pt.bg=c('gray30','gray80'))
abline(h=0)

#plot(PAdata$Growth.Rate~PAdata$N.from.grazing.perCperday)
#abline(a = 0,b = 1)

rainbowplot(PAdata[PAdata$Light>0,],'N.deficit','Unused N (gN per gN per day)','Light Level (grouped by strain)',spacer,strains,0,'topleft',logchoice='n',ylimits=c(-1,12))

Generally, numbers are positive, indicating that (so long as we assume 100% conversion efficiency) Ochromonas cells can meet their N demand from phagotrophy.

Nitrogen budget for the total culture

NO3.molar.media <- 0.000882  # concentration of N in Nitrate in K media
# https://ncma.bigelow.org/PDF%20Files/DOC-050.000%20NCMA%20Website%20algal%20medium%20K%202020.pdf


Approx.Och.Conc <- 200000 # approx max population size of Ochromonas per mL while still in exponential growth

PAdata$TotalNinOch.pgNpermL <- PAdata$N.perOch*Approx.Och.Conc
PAdata$TotalNinOch.molperL <- PAdata$TotalNinOch.pgNpermL * 1000 / (10^12) / 14.01 # 1000 mL per L; 10^12 pg per g; 14.01 g per mol N


bargraph.CI(Strain,TotalNinOch.molperL*1000,group=Bact.short, data=PAdata[PAdata$Strain!=2616,],las=1,legend=F,xlab='Ochromonas Strain',ylab='mmoles per L N in Ochromonas',leg.lab=c('High Bacteria','Low Bacteria'),col=mega.col.strains,err.width = .05,ylim=c(0,NO3.molar.media*1000))
legend(x='topright',inset=0.0,legend=c('High Bacteria','Low Bacteria'),bty='n',pt.cex=2,pch=22,pt.bg=c('gray30','gray80'))
abline(h=NO3.molar.media*1000,lty=2)
text(0,NO3.molar.media*950,'Nitrate added to media (mmoles per L)',pos=4)

We can see that the cultures used 10% or less of the total N available to them in the culture, suggesting that cultures were not N-limited at the time that physiological data were being collected.

LOCATION OF STRAIN ORIGINS

newmap <- getMap(resolution="low")
par(mar=c(0,0,0,0))
plot(newmap,xlim=c(-180,180),ylim=c(-90,90),col='white')
text(metadata$TextLon-5,metadata$TextLat,metadata$Strain,pos=4)
arrows(metadata$Lon,metadata$Lat,metadata$TextLon,metadata$TextLat,length=0,lwd=1.5)
points(metadata$Lon,metadata$Lat,pch=8,cex=1.5,col='red')

Supplemental figure on bacterial abundance

# Some data cleanup
bact.counts$Bact.short <- substr(bact.counts$Rice,1,1) # Uses the first character of the bacterial treatment as the identifier
bact.counts$Strain <- factor(bact.counts$OchStrain,levels=strains) # Orders the strains for plotting according to growth differential

par(mar=c(4,5,1,1))
bargraph.CI(Strain,CFUpmL,group=Bact.short,data=bact.counts[bact.counts$OchStrain!=2616&bact.counts$OldData!='Y',],log='y',ylim=c(2e5,53000000),las=1,xlab='Ochromonas Strain',ylab='', col=mega.col.strains,err.width=.05)
mtext('Bacterial CFU per mL', side = 2, line = 4)
legend(x='topright',inset=.0,legend=c('High Bacteria','Low Bacteria'),pch=22,pt.cex=2.3,horiz=T,bty='n',pt.bg=c('black','gray80'))
x.coords <- c(2,5,8,11,14,17,20,23)

text(x.coords,c(2.8e7,3.5e7,1.4e7,4.2e7,4e7,2.5e7,2e7,2.8e7),c('***','***','***','***','***','**','*','***'),cex=1)

Significance testing

for(i in 1:length(strains)){
  print(strains[i])
  print(t.test(bact.counts[bact.counts$OchStrain==strains[i]&bact.counts$Bact.short=='X',]$CFUpmL,bact.counts[bact.counts$OchStrain==strains[i]&bact.counts$Bact.short=='R',]$CFUpmL)$p.value*8)
}

## [1] 584
## [1] 3.039584e-06
## [1] 590
## [1] 1.734052e-10
## [1] 2951
## [1] 6.734511e-06
## [1] 1393
## [1] 9.339223e-07
## [1] 1392
## [1] 5.446116e-10
## [1] 1148
## [1] 0.004671076
## [1] 1150
## [1] 0.0368448
## [1] 1391
## [1] 8.732771e-08

Supplemental figure on ingestion rates

# Set the physical amount of space between each strain's plot
spacer <- 80

# Use the 'rainbowplot' function defined above to make the plot. Syntax: 
#   First argument = dataset
#   Second argument = y-axis variable
#   Third argument = y-axis label
#   Fourth argument = x-axis label
#   Fifth argument = amount of space between plots
#   Sixth argument = order of the strains
#   Seventh argument = horizontal line position, or NA if none needed
#   Eigth argument = location for the legend
#   Ninth argument = log scaling of y axis, enter 'y' or 'n'
#   Tenth argument = y axis limits (optional)
#   Eleventh argument = expansion factor for x-axis label (optional)
#   Twelfth argument = 'n' if omitting top labels (optional)
#   xtickomit == omits x tick labels (optional)
rainbowplot(PAdata,'C.from.grazing.perCperday','','Light Level (grouped by strain)',spacer,strains,0,'topleft','y')

# y-axis label
mtext('Phagotrophy (g C per g C per d)',side=2,line=3.3)

Phylogeny explorations

How to conduct a Mantel test: https://stats.oarc.ucla.edu/r/faq/how-can-i-perform-a-mantel-test-in-r/

# excise 2616 data from matrix of base pair substitutions
gendist <- gendist[gendist$X!=2616,]
gendist <- gendist[,1:9]


# create a matrix of physiology summary data to perform correlational analysis on
strainsumm <- as.data.frame(gendist[,1])
colnames(strainsumm) <- 'Strain'

growth.summ.C$Strain3 <- factor(growth.summ.C$Strain,levels=strainsumm$Strain)

# insert max data
strainsumm$growth.max <- tapply(growth.summ.C$Growth.Rate,growth.summ.C$Strain3,FUN= max,na.rm=T)
strainsumm$chl.max <- tapply(growth.summ.C$mean.chl,growth.summ.C$Strain3,FUN= max,na.rm=T)
strainsumm$att.max <- tapply(growth.summ.C$mean.att,growth.summ.C$Strain3,FUN= max,na.rm=T)
strainsumm$phot.max <- tapply(growth.summ.C$mean.phot,growth.summ.C$Strain3,FUN= max,na.rm=T)
strainsumm$graz.max <- tapply(growth.summ.C$mean.graz,growth.summ.C$Strain3,FUN= max,na.rm=T)

# insert data ranges
strainsumm$growth.range <- tapply(growth.summ.C$Growth.Rate,growth.summ.C$Strain3,FUN= max,na.rm=T)-tapply(growth.summ.C$Growth.Rate,growth.summ.C$Strain3,FUN= min,na.rm=T)
strainsumm$chl.range <- tapply(growth.summ.C$mean.chl,growth.summ.C$Strain3,FUN= max,na.rm=T)-tapply(growth.summ.C$mean.chl,growth.summ.C$Strain3,FUN= min,na.rm=T)
strainsumm$att.range <- tapply(growth.summ.C$mean.att,growth.summ.C$Strain3,FUN= max,na.rm=T)-tapply(growth.summ.C$mean.att,growth.summ.C$Strain3,FUN= min,na.rm=T)
strainsumm$phot.range <- tapply(growth.summ.C$mean.phot,growth.summ.C$Strain3,FUN= max,na.rm=T)-tapply(growth.summ.C$mean.phot,growth.summ.C$Strain3,FUN= min,na.rm=T)
strainsumm$graz.range <- tapply(growth.summ.C$mean.graz,growth.summ.C$Strain3,FUN= max,na.rm=T)-tapply(growth.summ.C$mean.graz,growth.summ.C$Strain3,FUN= min,na.rm=T)

# compute the distance matrix for physiology
# All physiology
physiodist <- dist(prcomp(strainsumm[,2:dim(strainsumm)[2]])$x,diag=T,upper=T)
mantel.test(gendist[,2:9],as.matrix(physiodist),graph=F,las=1)

## $z.stat
## [1] 29209.13
## 
## $p
## [1] 0.201
## 
## $alternative
## [1] "two.sided"

# Maxima
physiodist <- dist(prcomp(strainsumm[,2:6])$x,diag=T,upper=T)
mantel.test(gendist[,2:9],as.matrix(physiodist),graph=F,las=1)

## $z.stat
## [1] 22137.37
## 
## $p
## [1] 0.181
## 
## $alternative
## [1] "two.sided"

# Range
physiodist <- dist(prcomp(strainsumm[,7:11])$x,diag=T,upper=T)
mantel.test(gendist[,2:9],as.matrix(physiodist),graph=F,las=1)

## $z.stat
## [1] 18655.04
## 
## $p
## [1] 0.239
## 
## $alternative
## [1] "two.sided"

#mantel1 <- mantel.test(gendist[,2:9],as.matrix(physiodist))

Results suggest no strong correlation between genetic distance and physiological differences.

Import the tree constructed by Matt Johnson

PhyloTree <- read.tree('/Users/hollyvm/GoogleSync/StudentWork/GinaBarbaglia/PhyloTree_rednames.nwk')
plot(PhyloTree)

require(phytools) # Package with phylosig function for significance testing of phylogenetic signal

## Loading required package: phytools

## Loading required package: maps

## 
## Attaching package: 'maps'

## The following object is masked from 'package:plyr':
## 
##     ozone

print('Chlorophyll')

## [1] "Chlorophyll"

## extract characters of interest
strain.chl.max <- setNames(strainsumm$chl.max,strainsumm$Strain)
## compute phylogenetic signal K
K.chl.max<-phylosig(PhyloTree,strain.chl.max,
    test=TRUE)

## [1] "some species in tree are missing from x , dropping missing taxa from the tree"

print(K.chl.max)

## 
## Phylogenetic signal K : 1.05102e-05 
## P-value (based on 1000 randomizations) : 0.182

strain.chl.max <- setNames(strainsumm$chl.range,strainsumm$Strain)
K.chl.max<-phylosig(PhyloTree,strain.chl.max,
    test=TRUE)

## [1] "some species in tree are missing from x , dropping missing taxa from the tree"

print(K.chl.max)

## 
## Phylogenetic signal K : 8.41103e-06 
## P-value (based on 1000 randomizations) : 0.153

print('Growth')

## [1] "Growth"

strain.chl.max <- setNames(strainsumm$growth.max,strainsumm$Strain)
K.chl.max<-phylosig(PhyloTree,strain.chl.max,
    test=TRUE)

## [1] "some species in tree are missing from x , dropping missing taxa from the tree"

print(K.chl.max)

## 
## Phylogenetic signal K : 5.01576e-06 
## P-value (based on 1000 randomizations) : 0.287

strain.chl.max <- setNames(strainsumm$growth.range,strainsumm$Strain)
K.chl.max<-phylosig(PhyloTree,strain.chl.max,
    test=TRUE)

## [1] "some species in tree are missing from x , dropping missing taxa from the tree"

print(K.chl.max)

## 
## Phylogenetic signal K : 8.88546e-06 
## P-value (based on 1000 randomizations) : 0.426

print('Attack')

## [1] "Attack"

strain.chl.max <- setNames(strainsumm$att.max,strainsumm$Strain)
K.chl.max<-phylosig(PhyloTree,strain.chl.max,
    test=TRUE)

## [1] "some species in tree are missing from x , dropping missing taxa from the tree"

print(K.chl.max)

## 
## Phylogenetic signal K : 0.00152984 
## P-value (based on 1000 randomizations) : 0.005

strain.chl.max <- setNames(strainsumm$att.range,strainsumm$Strain)
K.chl.max<-phylosig(PhyloTree,strain.chl.max,
    test=TRUE)

## [1] "some species in tree are missing from x , dropping missing taxa from the tree"

print(K.chl.max)

## 
## Phylogenetic signal K : 0.00302222 
## P-value (based on 1000 randomizations) : 0.013

print('Phot')

## [1] "Phot"

strain.chl.max <- setNames(strainsumm$phot.max,strainsumm$Strain)
K.chl.max<-phylosig(PhyloTree,strain.chl.max,
    test=TRUE)

## [1] "some species in tree are missing from x , dropping missing taxa from the tree"

print(K.chl.max)

## 
## Phylogenetic signal K : 1.43675e-06 
## P-value (based on 1000 randomizations) : 0.732

strain.chl.max <- setNames(strainsumm$phot.range,strainsumm$Strain)
K.chl.max<-phylosig(PhyloTree,strain.chl.max,
    test=TRUE)

## [1] "some species in tree are missing from x , dropping missing taxa from the tree"

print(K.chl.max)

## 
## Phylogenetic signal K : 1.59771e-06 
## P-value (based on 1000 randomizations) : 0.76

print('Graz')

## [1] "Graz"

strain.chl.max <- setNames(strainsumm$graz.max,strainsumm$Strain)
K.chl.max<-phylosig(PhyloTree,strain.chl.max,
    test=TRUE)

## [1] "some species in tree are missing from x , dropping missing taxa from the tree"

print(K.chl.max)

## 
## Phylogenetic signal K : 1.01605e-05 
## P-value (based on 1000 randomizations) : 0.405

strain.chl.max <- setNames(strainsumm$graz.range,strainsumm$Strain)
K.chl.max<-phylosig(PhyloTree,strain.chl.max,
    test=TRUE)

## [1] "some species in tree are missing from x , dropping missing taxa from the tree"

print(K.chl.max)

## 
## Phylogenetic signal K : 9.8714e-06 
## P-value (based on 1000 randomizations) : 0.384